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Abstract 



We propose a new renormalization scheme of the running couphng constant in 
general gauge theories using the Wilson loops. The renormalized coupling constant is 
obtained from the Creutz ratio in lattice simulations and the corresponding pertur- 
bative coefficient at the leading order. The latter can be calculated by adopting the 
zeta-function resummation techniques. We perform a benchmark test of our scheme 
in quenched QCD with the plaquette gauge action. The running of the coupling con- 
stant is determined by applying the step-scaling procedure. Using several methods 
to improve the statistical accuracy, we show that the running coupling constant can 
be determined in a wide range of energy scales with relatively small number of gauge 
configurations. 
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1. Introduction 



One of the key subjects upon which recent attention has been focused is the flavor 
dependence of SU{N) Yang-Mills theories. In particular, given a number of flavors Nf, the 
question is whether the theory has an (approximate) infrared fixed point. This question is 
triggered by efforts to construct an alternative mechanism of electroweak symmetry breaking, 
via assuming the existence of a new, strongly interacting sector beyond the electroweak 
scale.l^ The earliest model of this sort, the so-called technicolor,'^'® gives rise to a dynamical 
electroweak symmetry breaking by introducing a QCD-like sector scaled up to some TeV. 
While theoretically appealing, the simplest form of the technicolor model and its variants 
with QCD-like dynamics are ruled out or disfavored by electroweak precision measurements. 
However, the possibility of such mechanism with a non- QCD-like theorj^-EJ jg still open, 
and may provide observable signatures at the LHC. It is thus an important but challenging 
task to investigate the low-energy landscape of spontaneously broken, strongly interacting 
gauge theories.'I^ 

Among the theoretical tools at hand, the numerical approach to lattice gauge theories has 
made it possible to gain quantitative information about strong dynamics of gauge theories. 
The current understanding can be summarized as follows. A vector-like gauge theory, e.g. 
QCD, is known to exhibit confinement and dynamical chiral symmetry breaking for small 
number of massless fermions, Nf, in the fundamental representation of the gauge group. 
When A'^^ is just below the value N'^'^ at which the asymptotic freedom sets in, the theory 
is conformal (unbroken chiral symmetry, no confinement) in the infrared. Such a theory is 
believed to remain conformal down to some critical value A^^, where the coupling becomes 
strong enough and the transition to the confined chirally broken phase occurs. The range 
Nj < Nf < N'^^ is called the conformal window. 

It is thus essential to investigate strongly interacting gauge theories in a wide range of 
parameters, such as the number of colors, the number of flavors, and the fermion repre- 
sentations.'^ Several modern lattice studies in this research direction have recently been 
performed.'^"'^ In particular, the authors of Refs.'^''^^ performed calculation of the run- 
ning coupling constant using the Schrodinger functional scheme, and found evidence for an 
infrared fixed point in SU{?)) gauge theory with Nf = 12. However, it is important to study 
the running coupling constant in different renormalization schemes, in order to conclude 
that the fixed point is not an artifact due to a particular renormalization prescription but 
a physical one. For this purpose, we propose a new renormalization scheme which uses 
the Wilson loops as a key ingredient. Such a scheme is applicable to general gauge the- 
ories as long as the Wilson loops can be defined, and provides an efficient computational 
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method for lattice gauge theories. Specifically, a renormalized amplitude is defined as the 
ratio among the Wilson loops, namely the Creutz ratio, and its perturbative counterpart. 
The former can be evaluated non-perturbatively by Monte Carlo simulation, while the latter 
calculated analytically once the underlying theory is specified. By properly defining the 
non-perturbatively renormalized coupling constant, its scale dependence is extracted using 
the step-scaling procedure, i.e., from the volume dependence of the coupling.'^ 

Applying our scheme will provide not only an independent check on the extent of the 
conformal window, but also several computational advantages. The Creutz ratio can be 
obtained without 0(a) discretization errors, provided these errors are absent in the lattice 
action. This means our scheme is in principle free from any 0(a) systematic effect. Further- 
more, this scheme only involves simple gluonic observables, therefore does not introduce any 
particular kinematical setup which can deteriorate the discretization error or break chiral 
symmetry. Therefore it can be applied to simulations with dynamical fermions of any type, 
without restrictions on A'^^. For these features, this scheme may be an attractive alternative 
to the Schrodinger functional scheme or the twisted Polyakov loop scheme.'2^^2E^ 

Before performing calculations for the gauge theories with dynamical fermions, as a 
benchmark test, we apply this new scheme to the computation of the running coupling con- 
stant in quenched lattice QCD. The numerical calculation is performed using the plaquette 
gauge action with periodic boundary conditions. These boundary conditions are chosen for 
simplicity. Nevertheless, it results in effects of degenerate vacua known as the "toron" 
Our scheme can, however, be applied in principle to any choice of boundary conditions, 
such as twisted boundary conditions, which ensure no unwanted zero-mode contributions 
by inducing non-trivial background configurations. Adopting several methods to improve 
statistical accuracy, we can determine the running of the coupling constant in a wide range 
of energy scales with a relatively small number of gauge configurations. 

Another essential ingredient of our scheme is the perturbative calculation of the renormal- 
ization constant. This is performed analytically using zeta function resummation techniques, 
which prove to be quite convenient. First of all, zeta function techniques offer a natural 
method to study the analyticity (and regularity) properties of the perturbative counterpart 
of the Creutz ratio. In addition, some algebraic rearrangements of zeta functions, originally 
due to Chowla and Selberg,!^ allow us to recast the expressions in terms of analytic functions 
accompanied by some exponentially converging series, whose evaluation is almost trivial and 
requires little computer power. The zeta function methods we apply can be easily extended 
to any boundary conditions and to the case of the Polyakov lines.!^ 

This paper is organized as follows. In the next section we give the definition of the new 
scheme. The perturbative calculation is illustrated in Sec. 3. Section 4 is devoted to the 
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details of our numerical simulations, after brief introduction to the step-scaling procedure. 
Section 5 contains discussion on the numerical results and comparison with other results 
in the literature. Finally, Sec. 6 summarizes our conclusions. The paper contains two 
appendices where technical details and simulation parameters are reported. Preliminary 
results of this work have been presented in Ref.'^ 

§2. Wilson Loop Scheme 

In this section, we define a new renormalization scheme, the ' Wilson loop scheme' . Let 
us consider an amplitude A whose tree-level approximation is 

= kgl , (2-1) 

where is the bare coupling constant, and /c is a coefficient of proportionality that does 
not depend on and can be explicitly calculated for a given underlying theory. With a 
non-perturbatively calculated amplitude A^^ at the scale /i, the renormalization constant 
= A^ ^ {ii) / A^^"^^ relates the renormalized coupling constant, to the bare one, 

leading to the relation 

/(/i) = . (2-2) 

Although in particle physics an S'-matrix element, i.e. a scattering amplitude, is usually 
adopted as A, to define g{^) one can equivalently use any physical quantity that can be 
perturbatively expanded and is proportional to g^ at the tree-level. 
We define the Wilson loop scheme by taking the 'amplitude' to be 

Aw{R;Lo;go) = '^'g^^^ {W{R,T- Lo,To))y^j,.^^^^^^ , (2-3) 

where W{R, T; Lq, Tq) is the Wilson loop with the temporal and spatial sizes T and R, on a 
lattice of the physical size Lq x Tq. In this work, we take Tq to be the same as Lq, and drop 
it in the argument of the Wilson loop. The scale Lq will be identified as the renormalization 
scale later. On a finite lattice, W, and thus Aw, also depend on the lattice spacing a which 
is determined by the bare coupling g^. The dependence of Aw on a is removed by taking 
the continuum limit, a — t- 0. A pictorial definition of the Wilson loop is shown in Fig. [TJ 
Using lattice perturbation theory (see Fig. [2]), Aw can be shown to be proportional to gQ at 
the lowest order. Thus, once the value of k is calculated, relation (12-21) leads to, after taking 
the continuum limit, a prescription to obtain the renormalized coupling: 
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Fig. 1. Wilson loop defined on the latticized space-time box. Tq, Lq and T, R represent the size 
of the box and the Wilson loop in the temporal and spatial directions, respectively; a is the 
lattice spacing. 

90 

W{R, T;Lq,Tq, a; Qq) = 

9o 

Fig. 2. Schematic illustration of the perturbative expansion of the Wilson loop. 

In the above expression we have made exphcit that in the continuum hmit A; is a regular 
function of R/Lq only. This will be proved in the next section. The remaining factor on the 
right hand side of Eq. (12 •4[) can be evaluated on the lattice as the Creutz ratio, 

^(R + 1/2; Lo/a) ^ - In f m^+l.T + l;Lo/a) WjR^tLM 

\W{R+l,T;Lo/a)WiR,T + l;Lo/a) 

where T = T/a and R = R/a. The value of x is evaluated by a Monte Carlo simulation. 
The renormalized coupling constant in the Wilson loop scheme can be written as 

The quantity g"^ depends on three different scales, Lq, R, and a; by taking the ratio to Lq, 
we use r = {R + a/2) / Lq, a/Lo, and Lq as the independent parameters. Fixing r to a specific 
value means fixing the renormalization scheme. The ratio a/L^ specifies the discretization 
of the box, and can be removed by taking the continuum limit, a/L^ 0. After fixing 
the two dimensionless parameters r and a/Lo, gl becomes a function of single scale, Lq. In 
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our scheme, following the step-scaling procedure, Lq is identified as the scale at which the 
renormalized coupling is defined. 

Ref.'^ offers an alternative definition of the QCD running coupling related to Wilson 
loops and discusses its finite-size effects in Lq/R. Differently from us, their scheme defines 
the coupling constant in the T — )■ oo limit, in which the renormalized coupling is related to 
the quark force. 

There are several advantages in using the Wilson loop scheme. An evident one is that 
our scheme does not contain 0{a) systematic effects as long as they are absent in the lattice 
action. This is because the Creutz ratio is free from 0(a) discretization errors, due to the 
automatic 0{a) improvements of the heavy quark propagator after the redefinition of the 
mass and the wave function.'S^J This is in contrast to the case of the Schrodinger functional 
scheme, in which the boundary counter terms give rise to additional 0{a) systematic errors. 
Such a particular kinematical setup also breaks chiral symmetry. Furthermore, this scheme 
only involves simple gluonic observables and can be easily applied to the case with any type 
of dynamical fermions without restriction to the number of fiavors. 

§3. Computation of k 

One of the indispensable ingredients of the scheme presented in the previous section is 
the calculation of the coefficient k in Eq. f l2-ip . It can be generically split into two terms: 

k = ko + ki, (3-1) 

where ko represents the zero-mode contribution, while ki can be expressed as 



ki = -2 Cf 



ORdT 



E' 



sm J 



(2vr)4^ \ no n 



2 



(3-2) 



T=R 



where the summation is taken over integer values of Ui {i=0,. ■ ■ ,3) except for the case no = 
ni = 71-2 = n-s = (indicated by the prime in the sum), and = + n1 + nl + nl. The zero 
mode contribution depends on the boundary conditions. In the following, we will concentrate 
on the case of periodic boundary conditions. In this case, ko was initially calculated in Kef^ 
For SU(3) gauge group, ko is given by: 

2 „ / R'^* 



The scope of this section is to present a method to compute the quantity k. The method 
we develop will be illustrated for the case of periodic boundary conditions, but it can be 
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applied, with minor changes, to the case of twisted or mixed boundary conditions. As we 
have seen, the contribution from the zero mode, ko, can be separated from the rest and is 
obviously regular. Thus to compute (and to prove the regularity of) k, we only need to 
consider ki. Our starting point is the quantity: 



S{T/Lo,R/Lo)= J2 



sin 



no 

no=— oo 



X 



oo 27rJ? 
Z-^ Z-^ T)f -I- -I- -nf -I- nf Z-^ 



'ni + ni + ni + ni ^ nl + nj + nl 

ni,n.2 13=1 u 1 z o ni,n2 u 1 z 



(3-4) 



ki can be obtained from S(T/ Lq, R/Lq) via 



Although it is not possible to find a closed form for S{T/ Lq, R/ Lq) in terms of elementary 
functions, the use of zeta function resummation techniques and basic analytic continuation 
allows us to recast S{T/Lo,R/Lo) into the form of a practically computable quantity, and 
to prove the regularity of k through an explicit calculation. The computation is carried out 
in a few steps. The first is the evaluation of the sum over by using the Poisson summa- 
tion formula. Then, the summation over rii and n2, is written in terms of the Epstein zeta 
functions. After these steps, the expression of S{T/Lo, R/ Lq) becomes compact. However, 
without further rearrangements, it is of little practical use. To this aim, it is convenient 
to rewrite the Epstein zeta functions using the Chowla-Selberg formula that renders the 
zeta functions into the form of elementary analytic functions plus some rapidly converging 
series. The subsequent step is to analytically perform the integrals introduced when using 
the Poisson summation formula, and finally perform the remaining summations numerically. 
Although the above procedure may seem involved, the actual implementation is rather sim- 
ple. The method has also the bonus of providing a proof of the regularity of the Creutz 
ratio, as we will explicitly show in the following. 

The first step of our procedure is to employ the Poisson summation formula: 



00 «QQ 00 «C 

E/(^3) = -^/(o)+ / dtm+2j2 / 

r^3=l ^ „=i^0 



_ /(t) cos{27mt)dt . (3-5) 

n3=l "^0 n=l "^0 

A straightforward application of the above relation to the function S(T/Lq,R/Lq) gives: 

°o sin °° r°° 

SiT/Lo,R/Lo) = 2 y y / cos{2n{m + R/Lo)t)Ctis; no)dt , (3-6) 

^ no ^ Jo 

no=— oo m=— oo " 
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where we have used the standard definition of the generahzed Epstein zeta function: 

OO CO 

Ct{s; no)= J2 J2 K + nl + nl + 1') . (3-7) 

ni=— OO 712=— OO 

The parameter s is a regulator, introduced to perform the necessary analytical continuations. 
The limit s — )■ 1 will be taken at the end of the calculation. It is interesting that the function 
S(T/Lo, R/Lq) can be entirely written in terms of the integral function 

POD 

Z{Q)= 1 cos{2nm) Ct{s]no)dt . (3-8) 



Although compact, the result Eq. f l3-6p requires further manipulation. A useful way to 
handle these functions is to make use of the Chowla-Selberg formula. Refs.'3Sl'''3S} develop 
the appropriate formalism that allows us to express C,t{s] uq) as the sum of analytic functions 
plus a rapidly converging series: 



r(s-l). o n.n „^ 27r °° ' 



p,q=—oo 



X {nl + t') "-'-'^'^ K,^, (^27ry/n2 + t2^/7T?) . (3-9) 
The other tool is the following integral formula (see Ref.>33): 

/>oo 

/ COS (27rr2t) + nl) K,., {2nip' + q^Y^'{t^ + n^) ^2) dt (3-10) 

Jo 



The procedure is now straightforward and consists in using the relations fl3-9p and fl3-10p in 
Eq. f l3-6p . Some computations lead to 

Z{f2) = y^^^^-ZA2r{s - 1) cos (7r((l - + 1/2)) (^) K^s-i)-i/2 {2nnno) 

p,q= — oo 

X K^s-i)-i/2 [271710 {p' + q' + . (3-11) 



It can be easily checked that, in the above expression, the limit s — ?■ 1 can be taken safely 
giving 



OO / 



= _^e-2-l^l"° + I Yl (/ + g' + r2')"'/'e-2-"ov/^W?^. (3. 12) 

p,q=~OD 
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It is evident that the quantity Z{Q) is regular. We can now substitute Eq. (13-121) into 
Eq. f l3-4p . and rearrange it as 

5(T/Lo, i?/Lo) = ^ [5o(i?/Lo) + 5i(T/Lo, R/ L^)] , (3-13) 

where we have separated the no = contribution from the remaining part which is exponen- 
tially suppressed. This separation leads to the definition 



where 



So{R/L,) = A,{R/Lo) + A^iR/Lo) + A^{R/L,) , 
S,{T/Lo, R/Lo) = B,{T/Lo, R/L,) + B^iT/L,, R/ L,) + B:,{T/L,, R/Lq) 



m=— oo 

H-oo oo 



2\m+iR/Lo)\ ' 

1 



m=— oo p,q= 
+00 oo 



\{p^ + q^ + \m + {R/Lo)\^f^ 



1 



m=-oop=l (r + |m+ (i?/Lo)P) 



~ sin g-2^|m+(i?/Lo)|no 
no- 

oo 



- TrrioT/Lo ^ 2|m + (i?/Lo) 

=1 ' m=— oo ' \ / / 



B2{T/Lo,R/Lo) = 2y V V — rrj 



oo „;„ 27rT„ +00 oo -27rnoy^p^>^+(R7MT 



, ^ sm —riQ ' „ 



3-14) 
3-15) 



3-16) 
3-17) 
3-18) 
3-19) 
3-20) 
3-21) 



Due to the exponential suppression, the terms Bi, B2, and -B3, and thus Si{T/Lq, R/Lq), are 
clearly regular. Therefore, to prove the regularity of k, we only have to show that So{R/Lo) 
is also regular. To show that the terms (13- 161) . f l3-17p . and fl3-18p also lead to a regular 
expression for Sq{R/Lo) (and to compute them), it requires further manipulations. 
The first term, (I3-I6p . can be computed analytically: 



A^{R/Lo) = ^ 



-^-^P(R/Lo) + ^P{-R/Lo) 



(3-22) 



where iIj{x) is the Euler psi function.!^ The remaining two terms, A2 and A3, can be 
rearranged by performing first the summation over m, and then by using the Chowla-Selberg 



10 



formula. We leave the details in Appendix 13 and present the results here: 



A2{R/Lo) = 8^cos(2j7ri?/Lo)A'o \ 2]Ti^p'^ + g2j + (/2-independent terms) , (3-23) 

oo 

A3{R/Lo) = 8 ^ cos{2j7TR/Lo)Ko {2j'Kq) + (/^-independent terms) . (3-24) 

<?,i=i 

Written as above, it is a trivial matter to see that, after taking the derivative of Ai, A2, and 
^3 with respect to /?, Sq{R/Lq), and thus k, is nicely behaved due to the exponential fall-off 
of the Kelvin functions K^{z). 

The last step of our procedure consists in evaluating the above expressions. The numerical 
computation of the sums in Eqs. f l3-19p . fl3-20p . fl3-2ip . ( l3-23p . and f l3-24p . does not present any 
problem due to the exponential suppression. Differentiating with respect to R, substituting 
T/Lq = R/Lq, and combining the results according to Eqs. fl3-2p . f l3-3p . (13 •4p lead to the 
result for k. Figure |3] shows the dependence of the function k with respect to R/Lq. Table 
[3] provides some indicative values for k. 




Fig. 3. The figure shows the dependence of k on R/Lq. The solid Hne represents k, in the contin- 
uum limit, according to the analytical expressions given in the text. The squares are results 
obtained using numerical lattice calculation for the sum in Eq. (|3-2p . As can be seen from the 
figure, the (continuum) limit of k for Lq/u — t- 00 exists and is finite. Also, the convergence of 
the lattice computation to the continuum value is faster for larger values oi R/Lq. 
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I, 

K 




I, 

K 




n 1 

U. 1 1 O\jo 


n '^n 


u. luuuo 


0.10 


0.17360 


0.35 


0.15694 


0.15 


0.17336 


0.40 


0.13970 


0.20 


0.17259 


0.45 


0.10885 


0.25 


0.17058 


0.50 


0.05556 



Table I. Values for k from the continuum calculation. 

§4. Numerical Simulation 

In this section, we will describe the details of our numerical simulations. For later use, 
we define the coupling-squared, 

~9l ^ kgl , (4-1) 

where r = + a/2)/Lo. Note that we express as a function of /3, r, and a/L^ instead of 
Lq, r and L^ja. The above redefinition is chosen for convenience, since /3, r, and ajL^ are 
the actual input parameters for the simulations. 

4.1. Step Scaling 

We begin by briefly reviewing the step-scaling procedure (see Refs.'^S'ES-ES for details), 
that we use to evaluate the evolution of the running coupling in a wide range of the energy 
scale on the lattice. 

The first step is to fix a value for r, and find a set of parameters, Lo/a), which produce 
the same value of for several different choices of Lo/a: 

{(/3f\(Lo/a)S^)),(/3(^\(Lo/a)(^)),---}. (4-2) 

We achieve this by tuning the value of /3 in such a way that the physical volume Lq is fixed 
for different values of Lo/a. We denote this fixed physical volume for the starting point of 
the step-scaling procedure by Lq. 

The next step is to vary the physical volume from Lq to sLq, which gives the evolution 
of the running coupling from the energy scale Lq ^ to {sLq)~^, where s is the scaling factor. 
This step can be performed by changing the lattice size from (Lo/a)*-^-' to s(Lo/a)*^^\ leaving 
each value of 13^^^ unchanged. Values of g^^ calculated with these new parameter sets should 
be considered as the coupling at the energy scale {sLq)~^ up to discretization errors, and the 
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extrapolation to the continuum limit can be taken 

1 



.2 ' 1 



sLq sL, 



(4-3) 



where Z is the renormalization factor as defined below Eq.f l2-ip . The resultant value of the 
coupling, gj^, should be considered as the renormalized coupling at the energy scale {sLq)~^. 
This is the way to obtain a single discrete step of evolution of the running coupling with 
scaling factor s. 

Next, we find a new parameter set for {(3^'^\ (Lo/a)^^-*), which reproduces the value of 
glj{l/sLo) obtained in the previous step. Here, we chose the parameter set in such a way 
that the new lattice size {Lq/o)^'^^ is equal to the original one, (Lo/a)^^^ From here, we can 
repeat exactly the same procedure described so far: we calculate g^^ with the parameter set 
s{LQ/a)^^^). By iterating this procedure n times, we obtain the evolution of the running 
coupling from the energy scale to (s"Lo)~^. 

4.2. Simulation Parameters 

We use the standard Wilson plaquette gauge action defined on a four-dimensional Eu- 
clidean lattice with finite volume Lq. In this work, we adopt untwisted periodic boundary 
conditions. However, it is straightforward to use other boundary conditions {e.g., twisted 
boundary conditions) when necessary. Gauge configurations are generated by using the 
pseudo-heatbath algorithm with over-relaxation, mixed in the ratio of 1:5. In the remainder 
of this paper, we use the word "a sweep" to refer to the combination of one pseudo-heatbath 
update sweep followed by five over-relaxation sweeps. In order to eliminate the infiuence 
of autocorrelation, we either take large enough number of sweeps between measurements, 
or adopt the method of binning with a large enough size of bin to estimate the statistical 
error reliably. We perform the numerical simulations based on the step-scaling procedure 
explained in the previous section for a fixed value r = 0.3. (The reason for this choice will 
be given in the next subsection.) We set the scaling parameter s = 1.5 with five different 
starting lattice sizes being L^/a = 10, 12, 14, 16 and 18, which means lattice sizes after the 
scaling at each step are L^/a = 15, 18, 21, 24 and 27, respectively. We take = 0.2871 
(which corresponds to gl, = ^f^Q g-, — 1.728) as the starting value of the first step of the 
step-scaling procedure. Tunings of the values of /3 (namely, finding values of /3 which satisfy 
gl (/3,r = 0.3, ^) = 0.2871 for each ^ = 10, 12, 14, 16 and 18 in the first step) are carried 
out by interpolating the data obtained from simulations for different values of L^/a and /3 
shown in Fig. HI Each data point in the figure is calculated from 200 gauge configurations 
with 1000-sweep separation between configurations. Once we obtain values of /3 which re- 
produce gl {(3,r = 0.3, ^) = 0.2871 for ^ = 10, 12, 14, 16 and 18, we carry out simulation 
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Fig. 4. 3^ obtained from simulations for different values of Lo/a and /3. 

for s = 1.5 step-scaling, namely simulations for = 15, 18, 21, 24 and 27 with the values of 
f3 we tuned. These results are used to take the continuum limit, then the resultant value of 
g'^ becomes a starting value for the next step. We iterate this procedure seven times. The 
combination of Lo/a and /3 used for the simulations are shown in Table ElS 

4.3. Simulation Details 

There are several practical steps to calculate the quantity g'^{(],r, L^/a) from numerical 
simulations. Here we explain various technical details of our computations. 

We use the APE smearin^^ of link variables defined by the following equation; 



= Projsuis) 



(4-4) 



where n and c denote the smearing level and the smearing parameter, respectively. The 



*' The values of (3 in Table II have numerical ambiguities coming from the statistical errors of the data 
used for the interpolation. These ambiguities propagate to the error of the physical scale at each step which 
we are trying to fix. However, it turned out that the effect of that fluctuation to the result of simulation 
for s = 1.5 at each lattice size was negligibly small compared to the statistical error of the simulation itself. 
Thus, we ignore those errors and resultant fluctuation of the physical scale in the rest of our analysis. 
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Step 1 


Step 2 


Step 3 


Step 4 


Lo/a 




Lo/a 




Lo/a 


/3 


Lo/a 


/3 


15 


8.31 


15 


7.80 


15 


7.44 


15 


6.968 


18 


8.25 


18 


7.83 


18 


7.45 


18 


7.040 


21 


8.27 


21 


7.86 


21 


7.49 


21 


7.076 


24 


8.32 


24 


7.91 


24 


7.55 


24 


7.156 


27 


8.40 


27 


7.97 


27 


7.61 


27 


7.243 



Step 5 


Step 6 


Step 7 


Lo/a 




Lo/a 


/3 


Lo/a 


/3 


15 


6.571 


15 


6.207 


15 


5.907 


18 


6.656 


18 


6.303 


18 


6.000 


21 


6.734 


21 


6.377 


21 


6.087 


24 


6.797 


24 


6.463 


24 


6.170 


27 


6.871 


27 


6.546 


27 


6.229 



Table II. Parameter sets, Lo/a and /3, used for the simulation. 



smearing is done for links in all four directions. The result does not depend on the value of 
c significantly, and we take c = 2.3 in the present study. Here, we need to find the optimal 
values of r = '^"^"^^ and the smearing level n, by considering the following requirements. For 
better control of discretization error, it is preferable to choose a larger value of r. Meanwhile, 
for the purpose of reducing the statistical error, it is better to take a smaller value of r and 
higher number of n. Fig. [S]shows the smearing-level dependence of in the case of /3 = 8.25 
and Lo/a = 18 as an example. From this figure, we find the statistical error is notably 
reduced even at the smearing level one. In order to avoid over-smearing, n should be smaller 
than R/2. This condition leads to the lower bound, Lo/a > {4n + l)/(2r). We summarize 
the bound from this requirement in Table IIIII We observe (see Fig. for the example of the 
case Lo/a = 18 at /3 = 8.25) that the data of (i? + 1/2) = 1.5 and 2.5 in higher smearing 
level are not reliable because of over-smearing. By considering all the above requirements, 
we find that (r, n) = (0.3, 1) is the optimal choice. 

Once we fix the value of r (r = 0.3 in our current study), we need to estimate the value 
of gl, for non-integer R. We interpolate the value of using a quadratic function: 

f{R + 1/2) = Co + ci{R + 1/2) + C2{R + 1/2)^, (4-5) 

with interpolation ranges for each lattice size listed in Table IIVI An example is shown in 
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0.1 z_ 
1 



2 3 4 5 6 7 8 

^+1/2 

Fig. 5. The values of g'^ with statistical error for several values of {R+1/2) in the case of /3 = 8.25 
and Lo/a = 18. Data connected by solid, dotted and dashed lines denote the data with 0, 1 
and 2 smearing levels, respectively. 



Tab 



n 


r = 0.25 


r = 0.30 


r = 0.35 


n = 1 
n = 2 
n = 3 


Lo/a > 10 
Lo/a > 18 
Lo/a > 26 


Lo/a > 8.3 
Lo/a > 15 
Lo/a > 21.6 


Lo/a > 7.1 
Lo/a > 12.8 
Lo/a > 18.5 



e III. The lower bound on Lo/a to avoid over smearing. 



Fig. E] where {R + 1/2) 
Lo/a = 18. 



5.4 corresponds to the interpolation to r = 0.3 in the case of 



Lo/a 


^+1/2 


Rfnin 


Rmax 


Lo/a 


i? + l/2 




Rmax 


10 


3.0 


2 


4 


18 


5.4 


4 


6 


12 


3.6 


2 


4 


21 


6.3 


5 


7 


14 


4.2 


2 


5 


24 


7.2 


5 


7 


16 


4.8 


3 


5 


27 


8.1 


6 


8 



Table IV. Ranges used to interpolate the value of g^. The column "72 + 1/2" is the value that 
corresponds to r = 0.3. 



The last step of the calculation is to take the continuum limit of from data obtained 
for different combinations of (3 and Lo/a listed in each column of Table [ITl We show two 
example plots in Fig. [6l which are the continuum extrapolations for Step 1 and Step 7. 
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0.002 0.003 

(«/Lo)' 



1.75 

0.005 



4— point linear - 

5— point linear 
5-point quadratic - 



(«/Lo) 



Fig. 6. The continuum limit of g'^. The left and right panels show Steps 1 and 7, respectively. 



Since our Wilson loop scheme does not contain 0{a) systematic errors, we extrapolate to 
the continuum limit using a fit function linear in (a/Lo)^. Four data points {L^/a = 27, 24, 21 
and 18) are used for this extrapolation (shown as red lines in Fig. [6]), and the resultant value 
is adopted as the central value of in the continuum limit. We also take the continuum 
limit by using a fit function quadratic in {a/ LqY with five data points {L^/a = 27, 24, 21, 18 
and 15) (indicated by pink curves in Fig. [6]), and the difference between the central values of 
two fits are adopted as the systematic error coming from possible higher order discretization 
effects. In Fig. El we have also plotted extrapolation by a linear function with five points of 
data for comparison. In this figure, resultant values of the continuum limit obtained from 
different fit functions are plotted at {a/L^f' = 0. (For better visibility, we slightly displaced 
the data obtained from 5-point quadratic and 5-point linear extrapolations.) All the error 
bars shown in Fig. |H] are statistical only. 

4.4. Numerical Results 

We now show the results of our simulations which were performed using parameters in 
Table [III with procedures explained in the previous section. Details of parameter choice and 
numerical results are summarized in Appendix [Bl 

4.4.1. Running coupling 

In Fig. [71 we plot the resulting values of g"^ and their statistical errors for Lo/a = 18, 
21, 24 and 27 for Steps 1 ~ 7. The continuum limit was taken in the way explained in the 
previous section, and both the statistical and systematic errors are estimated. In Fig. [7] the 
values of g"^ in the continuum limit are shown with statistical and systematic errors added 
in quadrature. 

The running coupling constant g^, is extracted by dividing g"^ by k{r = 0.3) = 0.1661. 
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1 











0.0005 


0.001 


0.0015 




0.002 0.0025 


0.003 



(«/Lo) 

Fig. 7. Results of simulations and continuum limit of g'^ in Steps 1 ~ 7 (from bottom to top). 



The evolution of the running coupling constant is obtained by connecting the resultant values 
for Steps 1 ~ 7 by assigning appropriate scales to these steps. We plot the results in Fig. [HI 
We define the starting energy scale of Step 1 as 1/Lq, and the evolution of the running 
coupling constant is plotted as a function of energy in units of 1/ Lq. In this figure, errors are 
accumulated with the evolution of the running coupling appropriatelj-l^ in the same way as 
explained in Ref.!^ . For comparison, we also plot scheme-independent perturbative running 
couplings with one-loop and two-loop approximation as well as three-loop approximation in 
the MS scheme (from bottom to top). In the high energy region, where the perturbative 
computation is reliable, the Wilson loop scheme is consistent with the perturbation theory. 
The figure also shows that our simulation is reaching deep into the low energy region, in 
which perturbative calculation is no longer reliable. 



*' For the appropriate procedure of accumulating error, we need values of derivative of the step scaling 
function a(u). Here, we used the result of u'^ polynomial global fitting of <y{u) to obtain approximate values 
of a'{u). See the next sub-subsection for detail. 
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(I)' (!) (fl (I)' (!) (if f 1 wl„ 

Fig. 8. Evolution of the running coupling constant in Wilson loop scheme, g'^, obtained from the 
step-scaling procedure. Horizontal axis shows the energy scale in units of 1/Lq. Three curves, 
from bottom to top, show scheme-independent perturbative running couplings with one-loop 
and two-loop approximation, as well as that with three-loop approximation in the MS scheme. 



4.4.2. Beta function 

From the results of the simulation, we can also extract the non-perturbative /3 function 
by using the method explained in Ref.l^ . To this end, it is useful to define the step scaling 
function in the continuum limit a^^ 

a{u)=gl{sL), u^gl{L). (4-6) 

We list the simulation results for the step scaling function in Table El In the week coupling 
region, a{u) can be perturbatively expanded ad^ 

a{u) = u + SqU^ + siu^ H , (4-7) 

with the coefficients 

so = 26olns, (4-8) 

si = (26olns)2 + 26ilns, (4-9) 

where 6o and hi are the one-loop and the two-loop coefficients of the /3 function, respectively. 
In Fig. [9l we plot the data listed in Table |V] together with two-loop perturbative curve 
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u 


a[u) 


1.728 


1.97(6) 


1.97 


2.23(7) 


2.23 


2.68(5) 


2.68 


3.29(2) 


3.29 


4.33(14) 


4.33 


6.89(22) 


6.89 


11.5(3) 



Table V. Simulation results for the step scaling function a{u). The values in parentheses represent 
total errors in units of the last digits. 



o(w) 




1 2 3 4 5 6 

u 



Fig. 9. Polynomial fit (dashed and dotted curves) of the data for the step scaling function a{u). 
Fit ansatz and resultant values of fitting parameters for two curves are explained in the text. 
The two-loop perturabative curve (solid) is also plotted for comparison. 

for (j{u) (solid), as well as two curves which are the result of the following two kinds of 
polynomial fit. For the upper (dashed) curve, we fitted u + sqU^ + Siv? + S2u'^ to the 
data in the range of m < 5, and obtained S2 = 0.0019(3). For the middle (dotted) curve, 
we fitted u + squ^ + siv? + S2U^ + s^^u^ to all the data, and obtained S2 = 0.0033(6) and 
S3 = —0.00048(9). As is expected, in the week coupling region, the data is well explained 
by the two-loop perturbative result, and polynomial functions fit well to the data also. 
Meanwhile, the figure clearly shows that neither two-loop perturbative curve, nor simple 
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|3(a) 




a 

Fig. 10. Non-perturbative (3 function in the Wilson-loop scheme. One-loop (green dashed curve) 
and two-loop (blue dotted curve) perturbative (3 functions are also plotted for comparison. 



polynomial fits can explain the behavior of the data in larger values of u. This is nothing 
but an indication of the emergence of the non-perturbative effect. 

The formula to obtain non-perturbative /3 function from the step scaling function was 
given in Ref.l^ as 




(4-10) 



By applying this formula recursively, we obtained the discrete /3 function as shown in Fig. [TOl 
Here, while the values of the (t{u) are directly given by the results of the simulatioio, we 
need to use the result of the fitting to obtain the values of a'{u). To obtain the latter, we 
adopted the global fitting explained above. As we noted, it fails to fit to the data in larger 
values of m, however, as can be seen from the figure, we can expect that it gives numerically 
approximate values. Also, we used the value of perturbative two-loop approximation for the 
value of P{u) at the smallest u as an initial input value to use Eq. fl4-10p recursively. In the 
figure, we also plotted the one-loop and two-loop perturbative /3 function for comparison. 
From this result, we confirmed again that for smaller coupling region, the data consistently 
reproduce the perturbative picture, while the non-perturbative effect begins to emerge for 
larger coupling region. 

*^ Unlike was done in Ref.,'^ we do not have to carry out interpolation in u since we have data exactly 
on the value of u we need. 
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5. Discussion 



We have hitherto concentrated on how the couphng constant runs under the relative 
change of the scale without referring to its absolute value. The absolute scale can be esti- 
mated using the Sommer scale, ro, defined by 

r^F(r)\ =1.65. (5-1) 

Based on the phenomenological potential models, tq corresponds to about 0.5 fm. In this 
work, however, we use ro only to set the reference scale for comparison with other scheme. 
The formula relating /3 and ro is given in Ref.,!^ 

ln(a/ro) = -1.6805 - 1.7139(;g - 6) + 0.8155(/3 - 6)^ - 0.6667(/3 - 6)^ (5-2) 

which is valid in the region 5.7 < /3 < 6.57. For our lattices, /3 values in Step 7 in Table Hi] 
are in this range. In Table IVH we summarize r^ja calculated using Eq. ( l5-2p . and the 
corresponding tq/Lq. In Ref.,'^ the values of ro in the range 5.7 < P < 6.57 are determined 
with errors linearly increasing from 0.3% to 0.6%. These errors are mainly statistical. Since 
we are only aiming at a rough estimate of the absolute scale in this work, we assign the 
corresponding size of errors to the results in Table I VII By extrapolation to the continuum 
limit, we obtaiiH 

= 4.75 ± 0.26. (5-3) 

Lo 

Here, the statistical error is negligibly small compared to the systematic error. 



Lo/a 


/3 


ro/a 




15 


5.907 


4.543(17) 


5.174(19) 


18 


6.000 


5.368(22) 


5.096(21) 


21 


6.087 


6.196(27) 


5.041(22) 


24 


6.170 


7.040(33) 


5.012(23) 


27 


6.229 


7.677(37) 


4.858(23) 



Table VI. The Sommer scale at each (5 of the Step 7, estimated using Eq. (j5-2p . 



We can now estimate the A scale in units of ro. Since it is obvious from Fig. [H] that 
g"^ is well approximated by two-loop perturbative running coupling at high energy region, 
it is reasonable to estimate the scale A by using the value of gljil/Lo) from the following 

*' Here we take the continuum limit and estimate the systematic error in the same way as we did in 
subsection 14.31 
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two- loop relation^ between LqA and (7^(l/Lo), 



bl 



V^o + Ms,(l/^o) 

where = ll/(47r)^ and 61 = 102/(4-7r)^ are the one-loop and the two-loop coefficients of 
the /3 function of quenched QCD. By substituting the value q^j^I/Lq) = 1.728, we find 

Lo/l^-^°°P ^ 0.0399. (5-5) 

Combining this result with the value of vq/Lq in Eq. fl5-3p . we obtain the value of yl^L°°'^ in 
units of tq as 

ro^WL°°'' = 0-190 ± 0.010. (5-6) 

We also estimated, in the similar way as above, the value of ro/l^~^°°P in the case of SF 
scheme by using the data reported in Ref.l^ , and found the following result: 

roylgp^°°P = 0.301 ± 0.025. (5-7) 

By fixing tq as a reference scale, we obtain the following ratio of /lgp^°°^ to 

^2— loop 

-fL. = 1.58 ±0.16. (5-8) 

As a consistency check, we have also extracted the ground-state potential from our data 
at the largest physical volume via double exponential fits ^ then used the potential to 
estimate tq and Tc defined as 

r^F(r)\_ =0.65. (5-9) 

Our results on these quantities are well compatible with those obtained in Ref.!® 

We also estimated the value of A^p^""^ / A^^"^ without relying on the measurement of 
any low-energy physical quantity. This can be achieved by comparing results obtained from 
the same combinations of values of (3 and L^/a in the two schemes, with a fixed physical 
box size Lq which is much smaller than I/Aqqd- To this purpose, we have carried out 
simulations using exactly the same values of /3 and Lo/a as one of the data sets in Ref.'^ 
. In Table IVIH we list the values of the coupling constant in the SF scheme, g^p (which 
is denoted by in Ref.l^), and results of Agp^°°^ estimated from them by using Eq. (15 •4p . 



The values of gl^ and the corresponding A^^"^ are also listed in this table 



*) 



These values 



From the difference between two-loop and three-loop (3 function in the SF scheme, the error in the 
estimation to A due to the higher order effect is about 3%. This is reasonably small compared to other errors 
in this work. 
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1 .r^ I n 

ljqI a 


P 


ysF 




9w 




J 2 — loop ; J 2— loop 


10 


7.8538 


1.8776(93) 


0.0539(9) 


1.9427(68) 


0.06043(71) 


0.891(18) 


12 


7.9993 


1.8811(38) 


0.0542(4) 


1.8620(66) 


0.05232(65) 


1.036(15) 


14 


8.1380 


1.884(11) 


0.0545(11) 


1.8028(73) 


0.04667(68) 


1.168(29) 


16 


8.2500 


1.864(10) 


0.0525(10) 


1.7662(79) 


0.04331(72) 


1.213(30) 



Table VII. Values for and Lov1^~'°°p in the SF and the WL schemes for several sets of L^/a 
and /3. Results of ^gp are taken from Ref.^D . The values in parentheses represent statistical 
errors in units of the last digits. 



result from an interpolation procedure using the data points shown in Fig. |H The results 
for ylgp^°°''/yl^L°°^ at each Lo/a is also listed in the same table. A fit linear in a/L^ to these 
results gives the following continuum-limit estimation: 

y^2 — loop 

= 1.78 ± 0.07 (stat.) ± 0.04 (sys.) (5-10) 

Here the systematic error was estimated by the difference between values in the continuum 
limit with linear and quadratic extrapolations. 

A comparison between the results of Eqs. (15 •8p and fl5-10p shows the "universality" of 
the estimate of ylgp'°°''/yl^L°°'', namely, the two different estimates are consistent with each 
other 



From the theoretical point of view, this "universality" might be a trivial result since 
quenched QCD has only one scale in the theory. However, from the numerical point of view, 
it is a rather non-trivial consistency check since one estimation involves the measurement of 
a low-energy physical quantity while the other is completely from high-energy physics. 

We notice that the error of (7^ is large partly due to the fact that we used the SF as 
inputs to set the scale for each /3. Therefore, the discretization error could be a combination 
(or the difference ) of those from WL scheme and SF scheme. This error may be dominated 
by the error from the WL-scheme since the effective physical scale is the Wilson loop size R 
but not Lq so that one expects 0{{a/R)^) error. 



*^ We notice that the g/Lq dependence of 5^ is rather large partly due to the fact that we used the 
results of the SF-scheme study as inputs to set the scale for each /3. Therefore, the discretization error could 
be a combination (or the difference ) of those from WL scheme and SF scheme. This error may be dominated 
by the error from the WL scheme since the effective physical scale is the Wilson loop size R but not Lq so 
that one expects 0{{a/R)^) error. 

**' Three-loop perturbative calculation (though quite challenging to carry out it in the WL scheme, while 
already done in the SF scheme), together with the improvement of systematic and statistical errors, would 
enable us to perform a even more precise test of this "universality" by defining the A scale with three-loop 
/3 function. 
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§6. Summary 



We proposed a new scheme for the determination of the running couphng on the lattice. 
Our method is based on the measurement of the finite volume dependence of the Wilson loop. 
Unlike the SF scheme, our method does not have any 0{a) discretization error, therefore 
the systematic effect arising from the extrapolation to the continuum limit is expected to 
be quite small. We showed results of numerical study for the quenched QCD as a feasibility 
test of our scheme. These results confirmed that our method led to the step-scaling of the 
coupling which is consistent with the perturbative running coupling at high energy. We 
also showed that the coupling calculated by this newly proposed scheme deviates from the 
two-loop approximation below a certain energy scale. This deviation arises from the effects 
that are not captured by the two-loop approximation. We have confirmed that our scheme 
works well for the calculation of the running coupling with relatively small number of gauge 
configurations, demonstrating that the statistical error is under control by properly choosing 
the smearing level and r. We expect that this new method is also applicable to the calculation 
of the running coupling in other gauge theories, including the SU{N) gauge theory with a 
large number of dynamical fermions, which will be studied in our future work. 
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Appendix A 

Formulas ^TM), and 



In the following, we will briefly show how to obtain formulas (13-231) . and fl3-24p . Basically 
the method we use is a repetition of the generic technique we have adopted in Sec. [31 we 
first perform the sum over m, and then use the Chowla-Selberg formula, which rearranges 
the expression as a sum of some analytic function plus a series, suppressed by the presence 
of the Kelvin functions Kyix). For the first term one has: 



hoo oo 



rt=— 

oo 



1 



m=— oo p,q 
oo 



p,q=l 

in 



; {p^ + q^ + \m + R/Lo\^) 



r(s) 



+ Yh){p' + <f)~' ' E^""' cos(2j7ri?/Lo)ir,_i (2j7r^ 

j=i 



/p2 _j_ q2 



E 

p,q=l 



r(s) 



+ 4^cos(2j7ri?/Lo)i^o (2jn^/^TV' 



(A-1) 



In the above formula, the first term disappears upon derivation with respect to R and thus 
will not contribute to k. The second term is exponentially suppressed due to the presence 
of Kq{z) and, thus, regular. An analogous procedure applies to A-s^R/Lq): 



+00 oo 

E E 



1 



m=— oo q 
oo 



-{q^ + \m + R/Lo\ 



E ^^rfi) + 4 E cos(2jvri?/Lo)i^o (2j7rg) . 



(A-2) 

q=l ' \2/ qj=l 

As before, the first term is independent of R and disappears, when differentiated with respect 
to R. The second term is exponentially suppressed and regular. 

Appendix B 

Summary of Simulation Parameters and Numerical Results 



Here, we summarize simulation parameters, and numerical results obtained from those 
simulations. In the tables below, A'conf, ^sweep and A^binsize respectively represent number of 
gauge configurations from which measurement were taken, number of sweeps between each 
configurations, and number of bin size when we estimate the statistical error. The values in 
parentheses in the column of gl, represent statistical errors in units of the last digits. We 
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also listed the values of g"^ in the continuum limit with the magnitude of total error at each 
step. 



Step 1 


Lo/a 




^conf 


-^swcep 


-^bin size 


9l 


15 


8.31 


200 


1000 


1 


0.2846(24) 


18 


8.25 


200 


1000 


1 


0.2999(24) 


21 


8.27 


400 


1000 


1 


0.3100(21) 


24 


8.32 


400 


1000 


1 


0.3119(27) 


27 


8.40 


200 


1000 


1 


0.3137(51) 


continuum limit : g"^ - 


= 0.328 ±0.010 



Step 2 


Lo/a 


/3 


Nconi 


-^sweep 


-^bin size 




15 


7.80 


200 


1000 


1 


0.3371(26) 


18 


7.83 


200 


1000 


1 


0.3465(29) 


21 


7.86 


200 


1000 


1 


0.3482(36) 


24 


7.91 


200 


1000 


1 


0.3585(48) 


27 


7.97 


200 


1000 


1 


0.3632(66) 


continuum limit : g^^ - 


= 0.371 ±0.012 



Step 3 


Lo/a 


P 


-^conf 


-'■ ' sweep 


-^bin size 




15 


7.44 


200 


1000 


1 


0.3831(31) 


18 


7.45 


400 


1000 


1 


0.4056(25) 


21 


7.49 


400 


1000 


1 


0.4125(33) 


24 


7.55 


400 


1000 


1 


0.4278(50) 


27 


7.61 


200 


1000 


1 


0.4249(92) 


continuum limit : g"^ - 


= 0.445 ± 0.009 
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Step 4 


Tjt\ / (7 




N r 
^ ' coni 


N 

^ ' sweep 


' bin size 


ilw 


15 


6.968 


10000 


1 


100 


0.4829(9) 


18 


7.040 


10000 


1 


100 


0.4959(12) 


21 


7.076 


10000 


1 


100 


0.5146(15) 


24 


7.156 


10000 


1 


100 


0.5204(20) 


27 


7.243 


10000 


1 


100 


0.5181(24) 




continuum limit : glj = 


0.547 ±0.004 



Step 5 






-^conf 


' sweep 


-^bin size 


9l 


15 


6.571 


10000 


1 


100 


0.6307(14) 


18 


6.656 


10000 


1 


100 


0.6489(19) 


21 


6.734 


10000 


1 


100 


0.6606(27) 


24 


6.797 


10000 


1 


100 


0.6794(33) 


27 


6.871 


10000 


1 


100 


0.6931(41) 




continuum limit : = 


0.719 ±0.024 



Step 6 


Lo/a 


/3 


-^conf 


' sweep 


-'^bin size 


9l 


15 


6.207 


10000 


1 


100 


0.978(3) 


18 


6.303 


10000 


1 


100 


1.016(5) 


21 


6.377 


10000 


1 


100 


1.069(6) 


24 


6.463 


10000 


1 


100 


1.075(7) 


27 


6.546 


10000 


1 


100 


1.074(9) 


continuum limit 


: 9l = 


1.144 ±0.037 



Step 7 


Lo/a 


P 


-^conf 


' sweep 


-^bin size 


9l 


15 


5.907 


10000 


1 


100 


1.811(5) 


18 


6.000 


10000 


1 


100 


1.833(7) 


21 


6.087 


10000 


1 


100 


1.846(12) 


24 


6.170 


10000 


1 


100 


1.861(13) 


27 


6.229 


10000 


1 


100 


1.917(27) 




continuum limit : g"^ = 


1.914 ±0.042 
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